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solution  in  the  coarse  grid  may  differ  widely  from  the  corresponding 
approximation  in  a  finer  grid,  if  the  frequency  defined  by  the 
inhomogeneous  problem  is  close  to  one  of  the  eigenfrequencies . 

This  may  be  detrimental  to  the  aim  of  the  coarse  grid  step;  namely 
to  remove  from  an  approximation  in  a  fine  grid  long  waves  which 
cause  divergence  of  fine  grid  iterations.  This  is  remedied  by 
eliminating  long  wave  contributions  to  the  residual  directly  in  the 
next  finer  grid.  The  success  of  this  procedure  still  depends  upon 
the  choice  of  the  subspace  from  which  the  long  wave  corrections  in 
the  finer  grid  are  taken.  Even  this  modified  method  fails  if  the 
frequency  pertaining  to  the  inhomogeneous  term  is  too  close  to  an 
eigenvalue,  unless  the  subspace  contains  the  pertinent  eigenfunction 

In  dealing  with  the  eigenvalue  problem  one  must  update  the 
eigenfunctions  and  eigenvalues  in  each  iteration  step  and  in 
addition  remove  long  wave  perturbations  pertaining  to  eigenfunctions 
with  frequencies  lower  than  that  currently  under  evaluation. 

Multiple  eigenvalues  and  eigenvalues  that  lie  close  together  are 
treated  simultaneously. 
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SECTION  I 
INTRODUCTION 


The  present  report  describes  some  examples  for  Brandt's 
multigrid  method  (Ref.  1).  The  primary  purpose  is  to  gain  some 
direct  working  experience.  We  have  not  used  ready-made  programs 
which,  presumably,  are  available  by  now,  and  we  have  even 
disregarded  some  of  the  finer  points  of  Ref.  1.  The  examples 
chosen  are  simple  but  not  entirely  trivial;  we  seek  solutions  to 
the  Helmholtz  equation  for  frequencies  which  exceed  the  lowest 
eigenf requency  of  the  problem.  Cases  of  this  kind  are  included 
in  Brandt's  work,  at  least  in  principle;  but  additional  considerations 
are  needed  if  the  frequency  is  above  that  pertaining  to  the  lowest 
eigenvalue  and  in  particular  if  the  frequency  is  close  to  one  of 
the  eigenvalues.  In  this  respect,  there  is  an  element  of  novelty 
in  the  present  discussions.  To  obtain  results  which  are  completely 
satisfactory  in  the  latter  case,  one  needs  the  eigenfunctions 
which  pertain  to  the  eigenvalues  in  the  vicinity  of  the  driving 
frequency.  A  modification  of  Brandt's  procedure  by  which  such 
eigenvalues  and  eigenfunctions  can  be  computed  is  therefore 
included . 

The  work  was  initiated  because  of  its  possible  usefulness 
for  the  flutter  problem.  One  must  realize,  however,  that  for 
this  problem  the  boundary  conditions  at  a  large  distance  are  of 
a  nature  which  precludes  the  occurrence  of  standing  waves.  The 
eigenvalue  problem  as  such  has  no  direct  bearing  on  the  flutter 
problem.  In  other  respects  Brandt's  multigrid  approach  is  likely 
to  be  very  useful  for  problems  of  this  kind,  especially  if  the 
frequencies  are  not  low. 
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SECTION  II 

DESCRIPTION  OF  THE  SAMPLE  PROBLEM 


The  region  for  which  the  numerical  solution  of  the  Helmholtz 
problem  is  carried  out  is  square  (Fig.  1) .  Along  the  boundary  of 
the  region,  Dirichlet  conditions  are  prescribed.  The  fine  grid 
is  obtained  by  dividing  each  side  of  the  square  into  12  intervals. 
The  grid  size  of  the  coarse  grid  is  four  times  that  of  the  fine 
grid.  This  choice  is  not  optimal;  Brandt  suggests  a  ratio  of  2:1 
for  the  grid  parameters.  Brandt  envisages  the  use  of  more  than 
two  grids,  while  our  examples  use  only  two  grids.  If  there  are 
more  than  two  grids,  then  one  must  decide  when  one  should  go  to 
a  finer  or  to  a  coarser  grid.  In  this  regard  the  present  example 
fails  to  cover  all  aspects  of  Brandt's  procedure. 

The  Helmholtz  equation  reads 


fxx  *  %  ‘  ° 

where  y  (the  square  of  the  frequency)  has  a  fixed  value.  The 
difference  form  of  this  equation  is 


h 


IH.  I 

✓  V 


(2) 


where  i  and  j  give  the  numbering  of  the  grid  in  the  x  and  y 
directions,  respectively,  and  h  is  the  distance  between  adjacent 
grid  points.  The  boundary  conditions  are 

(3) 

0  =0  J  0  =  0,  o  *  /l,  j  =  0,  j  - 

■  V 

In  the  fine  grid,  one  has  121  unknowns  (the  values  of  j  at  the 
inner  grid  points) .  In  the  coarse  grid,  one  has  a  system  of  four 
equations . 


Because  of  the  lack  of  an  inhomogeneous  term,  this  system 
has  only  the  trivial  solution  .  =  0  unless  y  is  an  eigenvalue. 
The  values  of  <|k  j  which  arise  in  the  course  of  the  present 
computations  therefore  give  the  errors  immediately.  The 
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solution  process  ought  to  reduce  the  errors  to  zero.  In  a  more 
realistic  problem  the  error  never  appears  separately;  one  encounters 
approximate  expressions,  -^j/  which  are  the  sum  of  the  exact 
solution  and  the  current  error.  The  presence  of  errors  manifests 
itself  in  terms  of  "residuals"  obtained  by  substituting  the 
current  approximation  into  the  difference  form  (Eq.  (2))  of  the 
differential  equation.  In  our  discussions  we  are  interested  only 
in  solutions  to  the  difference  equations,  the  question  of  truncation 
error  is  of  a  different  nature. 
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SECTION  III 


EIGENVALUES  AND  EIGENFUNCTIONS 


In  the  present  example  one  can  determine  the  eigenvalues 
and  eigenfunctions  for  the  difference  form  of  the  Laplace  operator 
in  a  closed  form.  Of  course  this  information  should  not  be  used 
in  the  numerical  work;  it  is  used  to  show  for  which  values  of 
difficulties  may  be  expected. 


The  (nonnormalized)  eigenfunctions  are  given  by 

fp  =  <tcyi{£>nso)  An/—  n>  =  /-  ■  ■  // .  »*.*/•■•// 

Ttfi,  '  /t  '  1  /?  J  j  j  '  j  j  J 

Substituting  this  into  Eq .  (2)  one  obtains 


(4) 
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h*-  * 
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/ k,'<C  i  */  m  t  i  ? 
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SWI,  =■  ^  '  '  '  # 


(5) 


The  eigenvalues  pertaining  to  Eq.  (1)  are 

(ft*  Ms*7  ~  (jt )  l  ^  j  V- 4V"  (6> 

The  two  expressions  approach  each  other  if  n  and  m  are  sufficiently 
small . 


For  the  coarse  grid  the  eigenfunctions  are  again  given  by 
Eq.  (4)  with  i  =  4  and  8,  and  k  =  4  and  8.  The  pertinent  eigen¬ 
values  are 


6  ^ 


*>*-  -  «2 


(7) 


The  eigenfunctions  for  the  coarse  grid  can  be  found  simply 
by  inspection 


->V  =./ 

-»t-  -  / 

~  Kj  ~  =  fa*  ~  ' 

yyt-  —  2. 

f%¥  m  fa*  ~/  /  9%S  fys  “  ' 

■*  Z  j 

»*t  -  /  ' 

?**•&*•'  J  ---/ 

•«.  -  L 
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The  first  four  values  of  |ih 


n  1  ,  m  =  1 


n  1  ,  m  2 
a  2 ,  m  1 


coarse  i|t'  it  1 


0.  129 

0  .  2  9 


fine  arid  1 


0.1.1630 

0.33610 


m  .! 


d i 1 ferent i al 
oqunt. ion 


0.1 3708 


0 . 34269 


n 


0  .  1 7  9 


0. 9 39 90 


0.94831 


SECTION  IV 

COMPUTATION  IF  \i  IS  NOT  AN  EIGENVALUE 


The  computation  proceeds  in  the  following  steps: 

1.  Introduction  of  a  starting  solution  for  the  fine  grid. 

2.  Fine  grid  iterations  (in  our  examples  carried  out  by 
the  alternating  direction  method) .  Terminate  the  procedure  if 
it  converges. 

3.  If  the  fine  grid  iterations  converge  slowly  or  diverge, 
then  one  computes  the  residuals  for  the  fine  grid  and  transfers 
them  to  the  coarse  grid.  (According  to  Brandt  it  is  most  economical 
if  one  simply  evaluates  the  fine  grid  residuals  for  the  points 

of  the  coarse  grid.) 

4.  Computation  of  corrections  in  the  coarse  grid,  either 
by  iteration  or  by  direct  solution  of  the  system  of  equations  for 
the  coarse  grid. 

5.  One  transfers  these  corrections  to  the  fine  grid  by 
means  of  an  interpolation  and  adds  these  corrections  to  the  fine 
grid  values  of  f . ,  obtained  at  the  end  of  step  2.  Return  to 

1  K 

step  2, 

We  make  the  following  observation:  Steps  3,  4  and  5  are 
carried  out  in  order  to  remove,  as  well  as  possible,  long  wave 
contributions  which  appear  in  the  results  of  the  fine  qrid 
iterations.  This  purpose  may  not  be  achieved  by  the  procedure 
just  described  if  |i  is  close  to  one  of  the  eigenvalues  of  the 
problem. 

This  difficulty  can  be  ascribed  to  the  fact  that  the 
eigenvalues  are  not  the  same  in  the  coarse  and  in  the  fine  grid. 

The  following  example  shows  that  the  corrections  obtained  in  the 
coarse  grid  may  be  quite  unrelated  to  the  corrections  needed  in 
the  fine  grid.  Assume  that  ii  coincides  with  a  coarse  grid 
eigenvalue.  The  corrections  in  the  coarse  grid  are  obtained 
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by  solving  an  inhomogeneous  system  of  simultaneous  linear 
equations.  If  u  is  a  coarse  grid  eigenvalue,  then  the  solution 
will  be  infinite  although  the  contribution  of  long  wave  errors 
in  the  fine  grid  is  definitely  finite. 

The  elimination  of  long  wave  contributions  can  actually 

be  accomplished  without  solving  the  problem  in  the  coarse  grid. 

One  notices  that  the  coarse  grid  corrections  depend  upon  a  number 

of  parameters  which  is  equal  to  the  number  of  grid  points  in  the 

coarse  grid.  By  an  interpolation  one  obtains  a  fine  grid  correction 

for  each  of  these  possible  coarse  grid  corrections.  The  elimination 

of  long  wave  contributions  in  the  fine  grid  can  therefore  be 

accomplished  by  adding  to  the  current  fine  grid  approximations  a 

linear  combination  of  the  interpolated  coarse  grid  corrections. 

The  coefficients  with  which  these  functions  are  multiplied  are 

determined  in  such  a  manner  that,  as  far  as  possible,  the  long 

waves  are  eliminated  from  the  fine  grid  approximation.  The 

criteria  by  which  these  coefficients  are  determined  must,  of 

course,  be  derived  from  the  fine  grid  residuals  (see  the  remark 

at  the  end  of  Section  II) .  Accordingly,  one  computes  in  advance 

the  residuals  in  the  fine  grid  belonging  to  the  interpolated 

coarse  grid  corrections.  In  the  present  case  there  are  four 

linearly  independent  residual  functions  of  this  kind.  Assume 

that  one  has  carried  out  a  number  of  alternating  direction  steps. 

One  then  possesses  a  fine  grid  approximation  and  pertaining  to 

it  a  fine  grid  residual.  From  this  residual,  and  subsequently 

from  the  approximation  for  ,  one  wants  to  remove  long  wave 

contributions.  Let  R(4>)  be  the  residual  for  the  approximation 

obtained  at  the  end  of  the  alternating  direction  iteration.  Let 

fch 

R(4>  .)  be  the  residual  for  the  i  coarse  grid  expression 

’coarse,  1 

after  interpolation  to  the  fine  grid.  In  the  present  example 
i  =  1 , . . . , 4 .  According  to  the  idea  just  described  above,  we  form 

*  “  Wetomi)  (8) 

c  ■> 

It  is  our  goal  to  determine  the  coefficients  so  that  long 
wave  contributions  are  (approximately)  removed  from  R. 
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This  task  could  be  done  perfectly  if  exact  long  wave 
eigenfunctions  were  known.  The  long  wave  eigenfunctions  are 
approximated  by  linear  combinations  of  the  interpolated  coarse 
qrid  expressions  <1  .  .  (If  one  would  compute  the  coarse 

grid  eigenfunctions  and  then  proceed  to  the  fine  grid  by 
interpolation,  then  one  would  obtain  just  such  expressions.) 
This  leads  to  the  following  procedure.  Let  the  scalar  product 
of  two  functions  u  and  v  defined  at  the  grid  points  be  given  by 


r  LC  V  ]  -  IL  & 
*•  J  - 


r 

*x-  yv 


Now  we  impose  the  requirement  that 

/  d  ,  £  7  =  o  ,  i  s  /  ■  •  *  (10) 

This  leads  to  the  system  of  equations 

//  i*  y  -  O  (11) 

where  M  is  a  matrix  with  components 

^ 

a.  y.^ 

!  is  a  vector  whose  n  component  is  given  by  a  ,  r  is  a  vector 
whose  mfcl1  component  is  given  by 

-  1131 

The  matrix  M  is  symmetric.  Let  L  be  the  difference  operator 
m,n 

defined  by  Eq.  (2) .  Then  one  can  define  an  eigenvalue  problem 
(L  ~  \  I)  ^  -  O 

The  eigenfunctions  4>^  are  orthogonal  to  each  other 

(14) 

The  residual  pertaining  to  4^  is  given  by 

X'fyJ  mLttJ  ■  <15’ 

Now  assume  that  <t>  ■  is  represented  in  the  form 

coarse , l 
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Then  one  has 


n 

nt/  >v 

j 


v 

Tceente)  >*v  ) 


- 

i 


it  kt  <&**■>  A k  WtJ]  A. \ j 

Finally,  because  of  the  orthogonality  relations  Eq.  (14) 

~  ^  (fcwb,  {fcvh.  \  </>*,] 


Hence 


4U 


ft. 


Oi/.  «M/ 


(16) 


The  system  Eq.  (11)  arises  by  forming  scalar  products  of 

the  expression  Eq.  (8)  with  the  functions  d)  .  In  a  first 

coarse, m 

version  we  had  introduced  the  condition  that  the  scalar  product 

of  the  residual  with  itself  should  be  a  minimum.  This  leads  to  a 

condition  similar  to  Eq.  (10) ,  except  that  now  the  weight  functions 

are  given  by  R(<i  „,)  ,  rather  than  <}>  itself.  If  y 

coarse, m  Tcoarse,m 

is  close  to  the  square  of  an  eigenfrequency ,  then  one  of  the 
eigenvalues  X  is  close  to  zero,  and  the  contribution  of  the 
pertinent  eigenfunction  is  nearly  suppressed  in  the  weight  function 
R(v*  ^  =  _)  .  This  is  undesirable  for  the  purpose  of  the  computa- 

tion  is  to  remove  the  contributions  of  these  eigenfunctions. 

According  to  these  considerations,  we  replace  steps  3,  4  and 
5  by  the  following  procedure.  In  preparing  for  the  computations, 
one  determines  the  elements  of  the  matrix  M  (a  4  by  4  matrix  in  the 
present  example).  At  the  end  of  step  2,  that  is,  after  one  has 
decided  that  the  convergence  is  too  slow  one  evaluates  the  residual 
R ( cf> )  pertaining  to  the  current  approximation,  forms  the 
components  of  the  vector  r,  Eq.  (13),  and  then  solves  the  system 
Eq.  (11)  for  the  vector  a. 

With  the  values  of  cx^  so  obtained,  one  forms  a  corrected 
function  <p 

tcs»ec/uf  '  t  +  f  %  tierce  <17) 
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and  returns  to  the  beginning  of  step  2.  Notice  that  the  system 
for  the  has  the  same  dimension  as  the  system  of  equations  which 
one  would  obtain  in  a  direct  treatment  of  the  coarse  mesh  (step  4 
of  the  procedure  outline  above) . 

th 

One  notices  that  if  y  is  close  to  the  square  of  the  j  eigen- 

frequency,  say  n,  then  the  contribution  of  the  eigenfunction  to 

the  expressions  R(<J>  )  and  consequently  to  the  elements  of 

coa  rs0  /fi 

the  matrix  M  is  small.  (Of  course  the  contribution  of  this 
eigenfunction  to  the  inhomogeneous  term  R(4>)  is  also  small.) 
Consequently,  the  effect  of  this  eigenfunction  on  the  determination 
of  the  coefficients  ou  becomes  small  and  the  goal  of  the  computa¬ 
tion,  namely  the  approximate  removal  of  the  eigenfunctions  4>  •  from 
the  current  approximation  4>  (see  Eq.  17)  is  not  achieved.  The 
interval  of  values  y  in  which  this  failure  occurs  depends  upon 

the  character  of  the  functions  4>  1  s  for  this  determines 

coarse, n 

how  strongly  eigenfunctions  other  than  the  4*^  (here  k  =  1,...,4) 

occur  in  the  A  's.  These  other  eigenfunctions  will  then 

coarse, n 

play  the  dominant  role  in  the  matrix  M  and  lead  to  faulty  values 

of  a.  The  width  of  the  interval  y  in  which  the  approach  fails 

depends  upon  the  interpolation  formula  by  which  one  proceeds  from 

the  coarse  to  the  fine  mesh.  Best  results  would  be  obtained  if 

the  functions  <p  span  the  same  subspace  as  the  eigenfunctions 

cosrsG f n 

<4n  (in  our  case  n  =  1,...,4). 

These  observations  are  borne  out  by  our  computations,  which 
have  been  carried  out  for  different  interpolation  routines 
from  the  coarse  to  the  fine  grid.  In  routine  1  linear  interpolation 
has  been  used.  In  routine  2,  we  have  used  a  third  degree  inter¬ 
polation  formula  for  both  the  x  and  y  direction.  (This  is  possible 
in  the  present  example  because  the  eigenfunctions  arise  from  a 
product  hypothesis  and  because,  with  the  boundary  points 
included,  the  coarse  grid  has  only  4  points  in  the  x,  and  4 
points  in  the  y  direction.)  In  routine  3,  exact  eigenfunctions 
are  used. 
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Figure  2  gives  a  survey  of  the  results.  The  horizontal 

1/2 

axis  shows  the  values  of  u  '  (one  might  just  as  well  use  y  as 
independent  variable) .  The  eigenfrequencies  for  which  we  expect 
difficulties  are  shown  by  short  heavy  lines.  The  vertical 
directions  show  the  number  of  iteration  steps  to  convergence. 

If  convergence  is  not  attained  by  five  iterations  but 
might  have  been  attained  by  more  iterations,  then  the  solid 
lines  are  drawn  up  to  the  upper  end  of  the  graph.  Divergence  is 
shown  by  dashed  lines.  As  expected,  there  is  an  interval  around 
the  eigenfrequencies  where  the  method  fails  to  converge.  This 
interval  becomes  smaller  as  one  uses  a  better  interpolation 
formula.  If  one  uses  the  exact  eigenfunctions  to  remove  long 
wave  perturbations  from  the  residual,  then  one  has  convergence 
even  very  close  to  the  eigenvalues.  In  this  case,  the  size  of 
the  interval  depends  only  on  the  precision  with  which  the 
computations  are  carried  out.  In  all  cases  the  convergence  is 
slower  for  larger  values  of  y .  No  convergence  is  obtained  if 
y  is  too  large. 
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( 


THIRD  DEGREE  INTERPOLATION 


CORRECT  EIGEN  FUNCTIONS 


Figure  2.  Number  of  Iterations  to  Convergence  for  the 

Inhomogeneous  Problem  with  Different  Functions 
for  the  Interpolation  from  the  Coarse  to  the 
Fine  Grid. 
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SECTION  V 

EIGENVALUE  PROBLEMS 


The  eigenvalue  problems  for  the  Helmholtz  equation  and  the 
Laplace  equation  are  the  same  except  for  a  shift  of  the  eigen¬ 
values  by  a  constant  amount.  The  following  discussion  is  there¬ 
fore  restricted  to  the  eigenvalue  problem  for  the  Laplace 
operator.  Let 


[L  (u,)]  -  1C, 

Then  we  consider 


+  u, 

Ljk- 


u,  / 

l-Hyk 


u, 

t-'ji 


(18) 


[Lfy)  4  *  0 

with  the  Dirichlet  boundary  conditions  formulated  above. 


(19) 


The  computations  start  with  approximations  for  a  limited 
number  of  eigenfunctions.  Such  approximations  can  be  obtained 
from  a  coarse  grid  formulation  with  a  subsequent  interpolation 
to  the  fine  grid.  In  our  example,  approximations  to  the  first 
four  eigenfunctions  have  been  determined.  If  one  wants  to  compute 
the  lowest  eigenfunction  only,  then  it  will  probably  suffice  if 
one  derives  from  the  coarse  grid  an  approximation  for  only  the 
lowest  eigenfunction. 


The  discussion  includes  cases  where  a  number  of  eigenvalues 
lie  close  together,  for  this  happens  frequently  in  multidimensional 
problems.  In  the  present  example,  for  instance,  one  has 
coincidence  of  the  second  and  third  eigenvalue.  The  coarse  grid 
approximation  will  show  for  which  eigenfunctions  this  is  the 
case . 


The  number  of  eigenfunctions  for  which  approximations  are 
introduced  is  restricted;  it  cannot  exceed  the  number  of  grid 
points  in  the  coarse  mesh,  but  one  might  use  even  a  smaller  number. 
The  subscripts  i  and  j  which  will  appear  subsequently  refer  to 
these  eigenfunctions.  Let  S  be  the  subset  of  subscripts  for  those 
eigenfunctions  which  are  treated  simultaneously^ because  the 
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corresponding  eigenvalues  lie  close  together  or  coincide.  The 
approximations  to  the  eigenfunctions  designated  by  the  subset  S 
are  updated  in  each  iteration  step,  those  not  pertaining  to  S 
remain  unchanged.  For  the  latter  one  always  uses  the  best 
approximations  available.  In  particular,  one  uses  the  exact 
eigenfunctions  (usually  obtained  in  previous  phases  of  the 
procedure)  if  they  are  available. 

We  shall  denote  eigenfunctions  and  their  approximations 
by  ,  if  i  f'  S  and  by  <J>j  if  j  t  S.  In  our  example  we  have  for 
the  first  eigenvalue  i  =  1,  j  =  2,  3,4;  for  the  second  and  third 
eigenvalues  i  =  2  and  3;  j  =  1  and  4;  and  for  the  fourth  eigen¬ 
value  i=4,  j  =  1,  2,  3. 

Let  a  superscript  p  denote  the  iteration  number.  One  iteration 

step  leads  from  t)’.  ^  to  ip .  ,  i  c  S.  The  iteration  starts 

1  ( 0 )  1 

with  approximations  tp  .  obtained  by  the  interpolation  from  the 
coarse  to  the  fine  grid.  We  describe  immediately  the  general 
case  where  S  contains  more  than  1  element.  Single  eigenvalues 
are  obtained  by  an  obvious  specialization. 

The  following  steps  are  carried  out. 

1.  Choose  the  subset  S,  provide  approximations  for  the 

*  j  ,  j  i  S,  and  starting  approximations  ,  i  e  S. 

2.  Update  the  values  of  and  form  an  intermediate  update 

to  the  approximations  ,  i  c  S. 

3.  Provide,  separately  for  each  i,  i  e  S,  a  second  update 
by  approximately  eliminating  long  wave  contributions  due  to  the 

•'  j  ' s ,  j  t  S. 

4.  Smooth  out  by  iteration,  separately  for  each  i  e  s 

the  twice  updated  functions  ip .  .  After  a  number  of  these  steps 
this  process  gives  the  final  approximation  .  If  one  has 

convergence,  then  one  updates  the  eigenvalues  and  terminates 

the  computation.  Otherwise  one  returns  to  step  2. 

We  discuss  these  steps  in  detail  and  provide  the  necessary 
equations.  In  step  2  we  set 


(20) 


f  -  Z  A'/' 

r  us  ^ 

i. 

with  coefficients  B^  so  far  unknown.  In  general,  such  functions 
will  not  satisfy  Eq.  (18)  , 

,  _  -  t- 

L  y>  ~  0 

even  if  one  considers  y  and  the  RVs  as  arbitrary.  To  obtain 
conditions  for  y  and  the  R4s  we  apply  to  the  last  expression  weight 
function  ,  i^e  S.  Then  one  obtains  the  eigenvalue  problem 

( M0>  =  0  (21) 

where  $  is  a  vector  with  components  8.,  i  e  S.  The  elements  of 

(1)  (2)  1 
M  and  M  are  given  by 


tl 


(22) 


(For  the  solution  of  this  problem,  one  will,  of  course,  change 
the  subscripts.)  The  solution  of  the  eigenvalue  problem  gives 
the  updated  eigenvalues  and  eigenvectors  and 

subsequently  (from  Eq.  (20))  updated  approximations  ^ . 

In  the  present  setting  the  matrices  and  M  ^  are 

symmetric  (see  the  discussion  in  Section  IV) .  Because  of  lack 
of  precision,  there  may  arise  deviations  from  symmetry  if  one 
evaluates  the  matrix  elements  separately.  It  is  then  desirable 
to  make  the  matrices  symmetric.  If  the  eigenvalues  y^p+^  are 
different,  then  one  obtains  eigenvectors  which  are  orthogonal 
to  each  other.  If  some  of  the  y^p+^  's  coincide,  then  one  can 
construct  vectors  which  are  orthogonal  to  each  other.  Critical 
are  cases  where  eigenvalues  y^p+^  are  different  but  lie  very 
close  together.  The  eigenvectors  obtained  in  such  a  case  will 
be  different  and  orthogonal,  but  small  errors  in  the  matrices  will 
lead  to  eigenvectors  which  are  quite  different  (although  taken 
together  they  will  always  span  the  same  or  approximately  the 
same  subspace  of  the  8  space) .  It  is  then  possible  that  the 
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eigenvectors  change  considerably  from  iteration  to  iteration 
depending  upon  small  changes  in  the  matrices.  This  will 
introduce  a  numerical  instability.  For  eigenvalues  which  nearly 
coincide  one  must  therefore  override  the  automatic  determination 
of  the  eigenvectors  and  define  eigenvectors  within  the  pertinent 
subspace  which  change  very  little  from  iteration  to  iteration. 


We  illustrate  by  example  another  disturbing  phenomenon 
which  may  occur.  We  consider  for  this  purpose  the  computation  of 
the  second  and  third  eigenfunctions.  In  this  case,  one  encounters 
2  by  2  matrices  which  have  nearly  diagonal  form.  Let  us  idealize 
them  as  diagonal  matrices. 


n 


to 


// 

O 


it 


*<r 


n 

/ 

o 


u> 


n 

a- 


The  determinant  vanishes  if 


(M 


(>' 


6  "  ** 


T  - ; 


A 


=  0 


Hence 


ft, 

ft 


o > 

SL 

u> 


or 


<T‘- 


iC- 

rt> 

Ll 


ft 


To  obtain  the  eigenvectors  one  substitutes  these  values  of  u. 
One  obtains  in  the  first  case 

1  o  o 

ct> 


0 

0 


r  \ 

A 

A 

l  ✓ 

Hence,  from  the  second  of  these  equations 


A * 0 '  A m ' 

With  inaccurate  numbers  the  zeros  which  appear  here  will  be 
small  numbers,  which  are  uncertain  because  of  the  lack  of  precision. 
If  one  inadvertently  determines  the  ratio  of  (3^  to  $2  from  the 
first  equation,  then  one  will  obtain  a  nonsensical  result,  while 
the  second  equation  will  give  the  correct  result.  Of  course. 
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a  good  eigenvalue  routine  will  guard  against  this  occurrence.  A 
problem  may  arise,  however,  if  one  uses  an  ad  hoc  program.  In 
the  present  example,  the  difficulty  is  easily  avoided  if  one 
determines  the  ratio  of  0^  to  02  from  the  sum  of  the  two  equations 
which  one  obtains  for  the  different  choices  of  y  . 

If  two  eigenvalues  coincide,  then  one  obtains  for  the 
determination  of  eigenvectors 

<0  a1fAl 


This  system  of  equations  is  satisfied  by  any  linear  combination 
of  0^  and  02* With  inaccurate  information  about  the  matrix  elements, 
one  will  obtain  specific  vectors  0^  and  02,  but  the  form  of  these 
vectors  will  depend  upon  the  inaccuracies  and  may  vary  from 
iteration  to  iteration.  To  guard  against  this  occurrence  we  have 
chosen  f  ,  1  f  o  l 


/*•  and  /'■l.J 

if  we  found  by  a  test  that  the  eigenvalues  are  very  close  together. 

In  the  following  steps  the  functions  i/m  with  approximate 
eigenvalues  y^p+1^  i  c  S,  are  treated  separately  for  each  i. 

In  step  3  the  long  wave  contributions  due  to  the  approximate 
eigenfunctions  <+> j ,  j  £  S  are  approximately  removed.  We  set 

3  it.  +  Z.  «0.d-  (23) 

Ti  71  jJS  ^ 

The  oij's  are  determined  by  the  requirement  that  this  expression 


be  orthogonal  to  the  functions  4>  ,  (j^  £  S) 

.  d) 

M  <  *  r  -  o 


One  then  obtains 


Where  a  is  a  vector  with  j  1  component  otj2»  i2  ^  S  anc3  the  matrix 
elements  are  given  by 

Vj  j  J.tSj  Jit*  (25) 


[%■„  hi 
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For  the  actual  computation  a  renumbering  of  the  subscripts 
will  be  carried  out. 

In  our  program  we  have  actually  proceeded  in  a  different 
manner  because  of  a  somewhat  uncritical  analogy  to  the  procedure 
for  the  inhomogeneous  problem.  We  believe,  however,  that  the 
method  descr_bed  above  is  preferable.  We  have  formed  the 
residual  pertaining  to  Eq .  (23) 

(L  '  <Ti  )*i 

J 

and  obtained  conditions  for  the  by  postulating  that  this 

expression  be  orthogonal  to  the  j  ^  £  S.  With  the  values  of 

i.  so  found  one  then  computes  ^  from  Eq.  (23)  . 

In  step  4,  one  eliminates  short  wave  perturbations  by  an 

iterative  procedure  (in  our  case  by  the  alternating  direction 

method).  This  is  done  separately  for  all  values  of  i,  i  >  S. 

If  these  iterations  converge,  then  one  terminates  the  procedure 

and  proceeds  to  another  set  S  of  eigenvalues.  Before  the 

termination,  one  may  update  the  eigenvalues  y  by  the  method  of 

step  2  for  a  last  time.  If  the  convergence  is  slow,  then  one 

goes  back  to  step  2  using  the  expression  obtained  by  the 

alternating  direction  iteration  as  starting  approximation 
(p+1) 

'  i 
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SKC’TJON  VI 

NUMKR1CA1.  KXPKRIKNCKS  IN  TMK  TRKATMK.NT  OK 
KIC.KNVA1.UK  PROlVIiKM 

Tho  proooiiuro  ilosoriboil  in  Soot  ion  V  has  boon  i  nip  1  onion  t  oil 
uni  nil  a  t  hi  ill  iloqroo  intorpolat  ion  formula  to  prooooil  from  tho 
ooarso  t  o  t  ht>  t  i  no  arid.  Wo  found  that  it  nmvt'nii’s  lather  wo  1  1 
With  1  inoar  i nt orpol at  i on ,  t ho  prooodnro  han  boon  triod  only  for 
tin'  n  initio  oiaonvaluos  i  1  ami  i  4.  Tho  oonvoraoiu'o  for 
i  1  is  sat  i  start  ory,  but  dubious  for  i  4.  Tho  quost  ion 

whi't  hoi  tho  prooodnro  rooommondod  in  tho  prooodinu  soot  ion 
for  st  i*p  4  wotilil  bo  bonofioial  has  not  boon  oxplorod. 
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